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Abstract. We study spatially partitioned embedded Runge-Kutta (SPERK) schemes for partial 
differential equations (PDEs), in which each of the component schemes is applied over a different part 
of the spatial domain. Such methods may be convenient for problems in which the smoothness of the 
solution or the magnitudes of the PDE coefficients vary strongly in space. We focus on embedded 
partitioned methods as they offer greater efficiency and avoid the order reduction that may occur 
in non-embedded schemes. We demonstrate that the lack of conservation in partitioned schemes 
can lead to non-physical effects and propose conservative additive schemes based on partitioning the 
fluxes rather than the ordinary differential equations. A variety of SPERK schemes are presented, 
including an embedded pair suitable for the time evolution of fifth-order weighted non-oscillatory 
(WENO) spatial discretizations. Numerical experiments are provided to support the theory. 
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1. Introduction. The method-of- lines is a popular discretization technique for 
the numerical solution of time-dependent partial differential equations. In it, a spatial 
discretization is applied to the PDE, yielding an initial value problem consisting of 
a large system of ordinary differential equations (DDEs). These are evolved by some 
time-stepping method, for example, a Runge-Kutta method. 

The choice of a suitable time-stepping scheme may depend on a variety of con- 
siderations. Typically, schemes are chosen to give good linear stability and accuracy, 
although in some applications it is also useful or necessary to preserve the mono- 
tonicity or other properties of the true PDE solution. Often, the key properties for 
determining a suitable time-stepping scheme vary locally in space. Examples of such 
properties include the grid spacing, the coefficients of the PDE, the smoothness of 
the solution and the geometry of the domain. As a consequence, it is possible that 
a scheme that is effective in one portion of the domain is unsuitable or inefficient in 
another. It is therefore natural to consider the development and analysis of spatially 
partitioned time-stepping methods, in which different step sizes or different methods 
are used over subdomains. 

A class of methods that often benefit from spatially partitioned time-stepping 
are PDE discretizations with grid adaptivity. In regions where the solution is nons- 
mooth or exhibits rapid variation, fine grids are needed; other regions may be more 
efficiently discretized with coarser meshes. If all components of an ODE system are 
evolved using some explicit time-stepping scheme and a single time step-size At, the 
evolution of all components will be restricted by the stiffest components of the system. 
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Improved efficiency is often possible by considering time-stepping methods that vary 
their time step-size according to the local mesh spacing. Time-stepping schemes of 
this type are called muUirate schemes. The first multirate schemes for one-dimensional 
conservation laws were developed by Osher and Sanders in [14]. Their approach car- 
ries out forward Euler time-stepping with step-sizes that vary locally. More recently 
higher-order methods have been considered; for example, Tang and Warnecke [22] 
develop second-order multirate schemes based on Runge-Kutta or Lax-Wendroff type 
schemes. Other recent work includes that by Constantinescu and Sandu [2] where 
a simple construction algorithm is given to form a second-order accurate multirate 
scheme based on an arbitrary strong-stability-preserving (SSP) scheme of order two or 
higher. Notably, their approach preserves a variety of monotonicity properties, such 
as positivity and the maximum principle. As explained by Hundsdorfer, Mozartova 
and Savcenco [7] multirate schemes for conservation laws are either locally inconsis- 
tent (e.g., [14]) or lack mass conservation (e.g., [22]). Fortunately, order reduction 
due to such local inconsistency may be less severe than expected, due to cancellation 
and damping effects; see [7] for details. 

Whereas multirate methods use different step sizes, in the present work we focus 
on using different methods on different spatial subdomains. Specifically, we inves- 
tigate spatially partitioned embedded Runge-Kutta (SPERK) schemes applied to 
one-dimcnsional conservation laws 

ut + fiu),=0. (1.1) 

Here the term embedded refers to methods having the same coefficient matrix A (but 
different weights 6), which avoids the unnecessary duplication of computations that 
can occur when combining two unrelated time-stepping schemes. Two classes of 
SPERK schemes are considered: equation-based partitioning and flux-based partition- 
ing. We shall find that equation-based partitioning maintains the overall accuracy of 
the schemes composing the embedded pair. This approach is, however, not conser- 
vative and can lead to wrong shock propagation speeds when applied to hyperbolic 
PDEs. Flux-based partitioning is conservative, but experiences a theoretical order 
reduction in the local consistency. In practice, however, we shall find that the overall 
accuracy of the combined scheme is unaffected. 

We will see that the methods we study can be viewed in the general framework 
of partitioned or additive Runge-Kutta methods. However, such methods are usu- 
ally applied in a way that applies different numerical treatment to different physical 
processes, whereas the emphasis here is on different numerical treatments for differ- 
ent spatial domains. These approaches are based on different motivations, and some 
approaches that work well for the former fail for the latter (see Example 2.3). 

As an example, consider a convection-diffusion problem in which the Reynolds 
number varies spatially, such that the system is dominated by convection on one 
subdomain and dominated by diffusion elsewhere. The third-order, three-stage SSP 
Runge-Kutta scheme of Shu and Osher (SSPRK(3,3)) [20] might be desired for the 
convection-dominated regions, while a second-order Runge-Kutta-Chebyshev method, 
e.g., [24], might be preferred in diffusion-dominated regions. Unfortunately, neither 
scheme is particularly attractive on its own since SSPRK(3,3) requires small time 
steps when applied to diffusive problems, while the Runge-Kutta-Chebyshev method 
is unstable when applied to convective problems. On the other hand, a combination 
that applies each time-stepping scheme where it is best suited may provide better 
linear stability than either scheme on its own. An example of this kind is explored in 
Section 4. 



Spatially partitioned embedded Runge— Kutta methods 



3 



There are other situations where spatiaUy partitioned time-stepping schemes show 
a strong potential. For example, consider the evolution of a conservation law involving 
both shocks and smooth regions. The preservation of monotonicity properties may be 
the most crucial property near the shock, whereas high-order accuracy may be of pri- 
mary interest in smooth regions. Use of a spatially partitioned time-stepping scheme 
opens the possibility of simultaneously obtaining good accuracy and monotonicity 
in such problems. An example of this kind is explored in Section 5, where SPERK 
schemes are applied to fifth-order weighted essentially non-oscillatory (WENO) spa- 
tial discretizations [19]. The motivation here is that fifth-order SSPRK methods are 
complicated by their use of the downwind-biased operator [20, 18, 21, 16] while mono- 
tonicity and the corresponding SSP property is likely only useful in the vicinity of 
non-smooth features. In these regions WENO discretizations are formally third-order 
accurate [19]. Thus, in our approach, a fifth-order linearly stable Runge-Kutta scheme 
is used in smooth regions while an embedded third-order SSP Runge-Kutta scheme 
is used near shocks or other discontinuities. 

The paper unfolds as follows. Section 2 introduces equation-based partitioning 
and examines its conservation properties. This is followed by the introduction and 
analysis of flux-based partitioning. Errors and positivity properties for both classes of 
schemes are considered in this section. A variety of generalizations are given in Sec- 
tion 3. Section 4 gives some examples of SPERK schemes, and considers applications 
of the methods to a spatially discretized advection-diffusion equation. Section 5 con- 
siders SPERK schemes in the context of WENO spatial discretizations and a nonlinear 
partitioning step. In our approach, the WENO weights are used to select between 
an SSP Runge-Kutta scheme, and a high-order scheme chosen for its linear stability 
properties. Finally, Section 6 concludes with a discussion of some other application 
areas and some of our current research directions. 

2. Spatially partitioned embedded Runge Kutta methods. We begin by 
introducing equation- and flux-based SPERK schemes. An analysis of the accuracy, 
conservation and positivity properties of such schemes is also provided in this section. 

2.1. Equation-based partitioning. Consider a system of N ordinary differen- 
tial equations 

U'{t)^G{U), (2.1) 

typically arising as the spatial discretization of a PDE where each component in the 
solution approximates, for example, point values of the PDE solution Ui sa u{xi) at 
discrete points Xi,l < i < N . To apply an s-stage Runge-Kutta method, we first 
compute the stage values 

s 

yO) =t/" + At^a,fcG(r«), j = l,...,s. (2.2a) 

fc=i 

A standard Runge-Kutta method would then advance by one step using the formula 
jjn+i ^ jjn _^ AtJ2'j=i bjGiY^i')). Instead, let a different set of weights be applied at 
each point Xi, 1 < i < N , choosing between coefficients b or 6. This results in 



S S 

X? E b,g^{Y^'^) + (1 - xD E bMy^'^) 



(2.2b) 
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where gi is the zth component of G and 

„ J 1, if weights b are used to compute u"~^^ 



u+i (2-2c) 



I 0, if weights b are used to compute 

We write the coefficients of this method A, b, b using the same tableau notation that 
is employed for embedded Runge-Kutta methods: 



6^ . (2.3) 
b^ 

We refer to the embedded method (2.2) as a spatially partitioned time-stepping 
method because the "mask" x selects which scheme to propagate at each point in 
space, at each time step. In some cases, x depends on the solution values, although 
we do not explicitly represent this for notational clarity. 

2.1.1. Connection to partitioned Runge Kutta methods. Methods of the 
form (2.2) form a special subclass of partitioned Runge-Kutta methods [5]. Generally, 
methods in a partitioned RK pair may have different coefficient matrices A as well 
as different weights b. In our embedded approach, methods have fewer degrees of 
freedom available for their design, but they possess the advantage of automatically 
satisfying the "extra" order conditions for partitioned RK methods. That is, if each of 
the component methods {A, b) and (A, b) is accurate to order p, the SPERK method 
is also accurate to order p in time (see also Proposition 3.1 below). 

Method (2.2) partitions (2.1) by equation; we refer to this type of partitioning 
as equation-based partitioning. Because the zth equation corresponds to grid node Xi, 
1 < i < N, our approach also gives a spatial partitioning of a semi-discretized PDE. It 
is worth noting, however, that equation-based partitioning is a very general technique 
that does not require any correspondence to grid locations. 

2.1.2. Conservation. Many important physical phenomena are modeled by 
conservation laws, which in one dimension have the form 

"* + /("). =0. (2.4) 

In the numerical solution of (2.4), one should use a conservative scheme in order to 
ensure that shocks propagate at the correct speed. Typically, (2.4) is semi-discretized 
using a flux-differencing method: 

where fi±i are numerical approximations to the flux at a^^-i-i , 1 < i < A^. Integrating 
with a Runge-Kutta method gives 

^6,/,+ .(r(^))-^fo,/,_.(rW) 

where the stage values Y'^^\1 < j < s, are defined in equation (2.2a) above. This 
method is conservative since corresponding fluxes cancel out (except at the bound- 
aries) if we sum over the spatial index z, 1 < i < A^. 
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Applying method (2.2) gives instead 
„ At 



Ax 



This method is not conservative, since the flux terms at x^_i will not cancel if x"-i 
Xi ■ This can lead to solutions in which discontinuities move at incorrect speeds. Here 
we give two examples. 

Example 2.1. Consider the advection equation Ut + Ux = discretized in space 
by upwind differencing: 



(2.6) 



where gi{U) = —-^ {ui — Ui_i). Let G — [gi{U)] and discretize in time with equation- 
based partitioning: 



y(2) ^ [/" + AtG(y(i)), 

<+i = K+At [x7 ■ 9^{y^'^) + (1 - xD • 9.{y^^^)) 



(2.7) 



The scheme (2.7) is a SPERK scheme and can be represented in the tableau form 
(2.3) as 












1 


1 







1 










1 



Consider the Cauchy problem on the real line with a step Junction initial condition: 

u{x, f = 0) = 



1 x < 0, 
X > 0. 



Suppose Xi = i, —oo < i < oo, t^ = n, n = 0, 1, . . . and for the spatial partitioning 
take 



1 i + 2n<l, 
i + 2n>l. 



Direct computation shows that the numerical solution is 

u- = i^ z<2n, 
]0 i>2n. 



This is not conservative since the discontinuity moves 2 grid cells per time step, 
whereas in the true solution it moves 1 grid cell per time step. Note that this ex- 
ample employs a scheme that is unstable for the given CFL number; the stability of 
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the solution results from the judicious choice of combination of methods, the initial 
condition, and the spatial partitioning. □ 
Example 2.2. Next we consider Burgers' equation, ut + (^u^)^ = 0, with the 
step function initial condition 

u{x, t = 0) ~ 

We discretize in space with a finite difference flux- differencing scheme using fifth- order 
WENO interpolation. In time, we discretize using the SPERK scheme in Table 5.1 
and equation-based partitioning. We use a time step-size of At = 0.6Aa;, correspond- 
ing to a CFL number of 1.2. We take 




X{x,t) = 



O.OK u(a;,i) < 1.99, 

1 otherwise. 



This simple choice of x ensures that the SSP method is used near the .shock while the 
non-SSP method is used elsewhere. It also (unfortunately) ensures that the jumps in 
X are near the shock, maximizing the effect of conservation errors. The true shock 
velocity is I, but the numerical shock velocity converges to approximately 0.925. If we 
replace x above by 1 — Xj .shock moves instead too rapidly. □ 

2.2. Flux-based partitioning. In (2.2) we applied one of the RK schemes to 
each equation in the ODE system. In order to obtain a conservative method, we 
partition by the fluxes /i+i , < i < N , rather than by the equations. 

Suppose we are given two Runge-Kutta methods with identical coefflcient matri- 
ces A, but difl'erent weights b and b. Similar to equation-based partitioning, we first 
compute stage values according to (2.2a), which when applied to the flux-differencing 
semi-discretization (2.5) is 



,0) 



At 



Vr' = n-i - ^ E«^-^ - , J = 1, ■ • • (2.8a) 

k=\ 

However, instead of varying the time-stepping method by equation, we vary the 
method by flux. An application of this partitioning to (2.5) yields 



- - (1 - x--i) E^./.-i(>^'-''')) , (2-8b) 

where we have partitioned using the characteristic functions corresponding to cell 
edges x^^i (rather than grid points Xi) 

{1, if weights b are used for fluxes at , i , 
^ . '+5' (2.8c) 

0, if weights b are used fluxes at 

Figure 2.1 compares a pseudo-code implementation of the equation-based (2.2) and 
flux-based partitioning (2.8) for a three-stage SPERK scheme. 

As we show next, the flux-based method corresponds to an additive Runge-Kutta 
method instead of a partitioned Runge-Kutta method. 
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2.2.1. Connection to additive Runge Kutta schemes. Consider again the 
flux-differencing formula (2.5). Let $(f/) denote the vector of fluxes, with components 
(t>i{U) = fi^i{U),0 < i < N. We can write the flux-differencing method (2.5) in 
vector form: 

U'{t) = -2^D^, (2.9) 

where D is the A'' x {N + 1) differencing matrix with 1 on the superdiagonal and —1 
on the main diagonal. 

We split the flux vector based on the vector x = [xi,X3,..., Xn+^]'^- 

U'{t) = -^D^U) = ^^D[diag{xMU) + (/- diag(x))<i>(f/)], 

" -^DdiagixMU) - diag(x))<i>(C/), 



= G^{U) + G^^^{U), (2.10) 

whereGx(C/) = -^D dis^g{x)HU), Gi-^{U) - -^i5(/ - diag(x))$([/) anddiag(x) 
is the {N + 1) X {N + 1) diagonal matrix which has x as its main diagonal. Flux-based 
partitioning applies different Runge-Kutta methods to and Gi-^. Because we 
have chosen embedded pairs with identical coefflcient matrices A, the stage values 
Y^^\ I < j < s, are computed according to (2.2a) in the standard fashion. Different 
weights for G^ and Gi_^ are applied, however, according to the formula 



^6,G;,(y(-'-)) + ^6,Gi_,(r(^) 



(2.11) 



This approach (i.e., using (2.11)) may be viewed as applying an additive Runge- 
Kutta method [9] to (2.10). We emphasize that the methods we consider have identical 
A matrices and therefore form a special embedded subclass of the additive RK methods. 
The effect of the splitting introduced in (2.10) is crucial for understanding the order of 
accuracy of flux-based SPERK schemes. Theoretically, we might expect a reduction 
of order by one, a result we establish in Theorem 2.4. In practice, however, we 
observe that the order of accuracy of SPERK schemes is the minimum order of the 
two schemes for both equation- and flux-based partitioning. Section 5 gives some 
numerical experiments illustrating this property. 

2.3. Accuracy. In general, additive and partitioned Runge-Kutta methods may 
exhibit lower order convergence rates than either of their component methods. In 
order for the method to be fully accurate, additional order conditions relating the 
coefficients of the different component methods must be satisfied. Even then, in 
our context of spatial splitting, such methods can exhibit order reduction when the 
coeflicient matrices A differ. 

Example 2.3. Flux-based partitioning based on a characteristic function x gives 
the ODE system (2.10). The classical Godunov operator splitting is a first-order 
method which may be applied to this system. Applying this splitting yields 

U* = C/"-l-AtGx(i7"), 



C/"+i = U* +AtGi-^{U*). 
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# Program 1: Equation -based SPERK implementation 

# First compute the standard stage values 
y(i) = U 

for i in ranged, N): 

yf = u, - dt/dx(a2i(yW-y(^)^)) 
for i in ranged, N) : 

= u, - dt/dx(a3i(y«-y(i>,) + a,^(jf^ -jf\)) 

# Advance in time from u to unew using: 
for i in range (1 , N) : 

unew, = u, - dt/dx( X, {bi{lf^ + 62(yf'-yf_\) + b-,(jf^ 

. d-X,)(&i(yW-y«) - S,(yf)-y(^J . 63(yf' -y'!),))) 

# Program 2: Flux-based SPERK implementation 

# First compute the standard stage values as in lines 3-7 above 

# Then advance in time from u to unew using: 
for i in ranged, N) : 

unewi = Mi - dt/dx^ 

- (l-X,+ i)(^'iyr' - feyf' - 'iayf') - (l-X,-i)(^>iyl^i - 'iay*^! - bsya)) 

Fig. 2.1. Pseudocode implementation of the equation-based and flux-based SPERK algorithms. 
Here we apply one step of an explicit three-stage method to the first-order upwinding spatial dis- 
cretization of the advection PDE ut -\- =0. 



This can he written as an additive Runge-Kutta method with coefficient arrays 





















1 


1 


, 













1 










1 



Unlike all other pairs considered in this paper, this pair uses different coefficients 
aij. We apply them to the upwind differencing semidiscretization (2.6) of the linear 
advection equation with At — Ax and constant initial data, — l,^oo < i < oo. 
We partition the domain into two regions by choosing an integer J with Xi^ i = 1 for 
i < J and Xi-\ — for i > J. After one step with this method, the computed solution 
is 

for i < J, 
for i — J, 
for i = J + 1, 
for i > J +1. 

Despite the constant initial data, the local error is 0(1) in the maximum norm. □ 
The difficulty encountered in the example above is typical of methods in which 
A ^ A. In the next two theorems, we establish the accuracy of equation- and flux- 
based SPERK methods. 

Theorem 2.1. Suppose that the partitioned Runge-Kutta method (2.2) with co- 
efficients {A,h,b) is applied to the semi- discretization (2.1). Then the fully discretized 
system has a local order of accuracy equal to min(p,p) where p (respectively, p) is the 
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order of accuracy of the full discretization obtained by using the Runge-Kutta method 
with coefficients (A,h) (respectively, {A,b)). 

Proof. The full discretizations are one-step methods of the form 

The exact solution satisfies 

u(a;„t„+i) - Q,{uH{tn)) + 0{MP+^), (2.12a) 

u{x,,tn+i) = Q^{uh{t,,)) + 0(AtP+i). (2.12b) 

where Uh{tn) — [M(a::i, t„), u(a;2, t„), . . . ,u{xf4,tn)Y' is a vector of true solution values 
at the grid nodes at time t„. Here we have assumed that Ax is given by some 
prescribed relationship in terms of At, so that the error can be characterized in terms 
of At only. 

We now determine the local truncation error of the discretization: starting from 
the exact solution at time t„, the solution computed by the partitioned method is 

= XlQ^{uh(tn)) + (1 - x7)Qr{uh{tn)). 

Applying (2.12), we find 



u 



n+l 



xr • «(x„ t„+i) + 0(Atf+i) + (1 - xD • u{x,,t,,+,) + 0(AtP+i), 



= u(x„t„+i) + o(Ar-(p^p)+i). 

□ 

In order to prove accuracy of the flux-based decomposition approach, we need to 
know how accurately a Runge-Kutta method approximates the fluxes. 

Lemma 2.2. Suppose we are given an initial value problem (2.1) for U G M", 
a Runge-Kutta method {A,b) of order p, and a smooth function V : M" — )• M". Let 
W{t) = J^V{U{s))ds, VF" = W{t„) and compute W^"+i using (2.2a) and 

S 

^"+1 ^ w" + AtJ2bjV{Y^^'>). (2.13) 

This scheme approximates W(t) to order p, i.e. W"~^^ — VF(t„+i) + 0{AtP'^^) . 
Proof. Suppose an 0{AtP) Runge-Kutta method is applied to the system 

U' = G{U), 
W' = V{U). 

Stage values are given by 

s s 

^U" + At ^ ajfcG(r('=)), Z(J') = M^" + At ^ ajkV{Y^''^),j 1, . . . , s. 

k=l k=l 

The Z''^\j — l,...,s, are not used and would not be computed in practice. We 
advance by one step via 

s 

^ jjn At^6jG(r(^)), 

J = l 

s 
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Recall that the Runge-Kutta method gives an error that is 0{AtP). Applying this 
fact to the second equation gives the desired result. □ 
From Lemma 2.2 we obtain 

Lemma 2.3. Let F : M" — >• M" be a smooth function, and let {A,b), (A, 6) be an 
embedded RK pair with order p, p. Given an initial value problem (2.1) in U , let Y'^^'^ 
denote the stage values in (2.2a); then 

s s 

bjV{Y^^^) = ^ hjV{Y^^^) + ©(At™"^^'^)). (2.14) 



Proof. Applying Lemma 2.2 to both schemes gives 

s 

W{t) + Ai ^ 6jy(y(^')) = W{t + At) + ©(AiP+i), 

W{t) + Ai ^ 6jy(y(^')) = W{t + At) + 0(AiP+i). 

j=l 

Combining these gives the stated result. □ 

One must take care in applying Lemma 2.3 to a PDE semi-discretization, since 
the error constant appearing in (2.14) might involve a factor like Ax""", r > 0, which 
would grow as the spatial grid is refined. A detailed analysis of the propagation of 
the spatial discretization error within a Runge-Kutta step is beyond the scope of this 
work; we refer the interested reader to [17, 1]. In the following theorem we simply 
assume that the error constant in (2.14) is bounded as Aa; 0. Note that this is an 
assumption on each of the component RK schemes separately. Theorem 2.4 indicates 
that if each of the component schemes gives a solution free from order reduction (in 
the sense just described), then the embedded pair loses at most one order of accuracy. 

Theorem 2.4. Suppose that the flux-based spatially partitioned Runge-Kutta 
method (2.8) with coefficients {A, b, b) is applied to the semi-discretization (2.5) with 
At = 0{Ax). Furthermore, suppose that the error constant appearing in (2.14) is 
bounded as Ax — > when one takes V = fi±^- Then the full discretization of (2.5) 
has order of accuracy equal to min(p, p) — 1 where p ( respectively, p) is the order of 
accuracy of the full discretization obtained by using the RK method with coefficients 
{A,b) (respectively, {A,b)). 

Proof. The full discretizations are 




^6,/,+ .(yW)-^5,/,_.(r(^)) , 



^6,/,+ .(y(^))-^S,/,_.(r(^))|. 



Using Lemma 2.3 with V{u) — fi±i{u) gives 

^&,/,;±i(F«) = 5]S,/,±.(y(^)) +0(At--(P^P)). (2.15) 
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n+1 



Flux-based partitioning gives 

J = l 3 = 1 

-xuiz b3h-i{y^'^) - (1 - x^ . ) E hh-h^Y^'' 




where we have used (2.15). The first two terms on the right are just the solution 
computed by the scheme with weights bj , which is accurate to order p by assumption. 
Thus, inserting the exact solution in the last equation above and using the assumption 
that At = ©(Ax) gives the desired result. □ 

Remark 2.1. More generally, for a semi- discretization of an evolution PDE that 
is order q in space, one takes (for explicit methods) At = ©(Ax"*) and a factor of 
Ax^^ appears in the spatial discretization. The net effect is the same in that it may 
reduce the accuracy by a factor of At. 

Remark 2.2. It might be tempting to try to apply the theorem to the general 
additive Runge-Kutta method 



A 



with 



A 



(2.16) 



by first making a larger additive Runge-Kutta method that does share abscissae and 
.stage weights, namely 



(2.17) 



c 


A 







c 


A 





c 





A 


with 


c 





A 




6^ 












IF 



and then concluding that Theorem 2.1 applied to the latter implies the results on the 
former. Indeed the theorem does apply to (2.17); however, the implication does not 
follow because (2.16) and (2.17) are not the same method. To see this, note that 
(2.16) computes stages combining function values from the two schemes (according to 
x) whereas (2.17) computes each stage value using only one scheme. 

2.4. Positivity. The nonlinear stability properties of SPERK schemes are of in- 
terest as well. In this section we consider the positivity of flux-based SPERK schemes 
which are comprised of two explicit SSP Runge-Kutta schemes. 

We begin by reviewing some results. Of particular relevance to our derivations 
is that SSP Runge-Kutta schemes may be re-written in optimal Shu-Osher form [18] 
via the transformation provided in [3]. For the ODE system U' ~ G{U) this yields 
schemes of the form 

/ i-i \ j-i 
rO) =il-J2^3k) U" + 5](«,fcr('=) + At/3,fcG(r('=))), (2.18a) 

V fc=i / fc=i 



jjn+1 ^ 



(2.18b) 



12 



Ketcheson, Macdonald and Ruuth 



where all ajk > 0, J2k=i ^jk ^ 1 ^^d j = 1, 2, . . . , s + 1. 

If both sets of coefficients ajk, 13 jk are nonnegative, and forward Euler is positivity- 
preserving, then it may be shown that the Runge-Kutta scheme preserves positivity 
under a suitable time step restriction [20, 4]: 

Lemma 2.5. // the forward Euler method is positivity-preserving under the re- 
striction < At < AtpE, then the Runge-Kutta method (2.18) with [3jk > zs 
positivity-preserving provided 



At < CAtpE, 

where C is the SSP coefficient 

C = min {cjk :l<fc<j<s + l} where Cjk 



otherwise. 



We now show the corresponding result for flux-based SPERK schemes. 

Theorem 2.6. Suppose that the flux- differencing semi-discretization (2.5) of the 
one- dimensional conservation law (2.4) satisfies the positivity properties 

- > 0,/or all 6ti < At^ (2.19) 

Wi + ^/i-i(^) > 0' V all St2 < At2 (2.20) 

for all Wi > 0,1 < i < N . Further, suppose that a flux-based SPERK scheme is applied 
to (2.5), and that the schemes composing the embedded pair have SSP coefficients C 
and C. Then the full discretization is positivity-preserving for time steps satisfying 

< At < min(C,C)At* (2.21) 



where 




if At I = oo, 
if At2 = oo, 
Au+At, Otherwise. 



Proof. We give a proof for the case where Ati and At2 are finite; the proof for 
the two remaining cases is straightforward and follows in a similar fashion. 

Let x"i_i and 1 — x^.i specify the (non-negative) weightings of the two schemes 
at the cell boundaries at time step n. The update for U is given by 
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where A = a|i 7p = At^+Ata ' ^^'^ ^"'^ ~ Att+At2 • denote the fluxes associated 
with stage j by f':'!?i = /j+ i (F'-"'-') where F^^^ is the approximation of Uh{tn) at the 

«± 2 2 

jth stage. The first term wih be non-negative if 

s 

7p< + A^-6jff, (2.23) 

is non-negative. Transforming (2.23) to optimal Shu-Osher form yields 

7p I 1 - ^a,+ij I < +7p I ^as+i,jyY' - ^/3,+ij/^\ 

We observe that forward Euler is positivity-preserving when applied to the semi- 
discretization (2.5) provided At < At*. An application of Lemma 2.5 gives y*--'^ > 
0, 1 < j < s, provided the step-size restriction (2.21) is satisfied and J/" > 0, where 
vector inequalities are to be interpreted as being taken componentwise. Combining 
this result with At < CAt* and hypothesis (2.19) yields that the first term in (2.22) 
is non-negative. The proof follows by applying a similar analysis to the second, third, 
and fourth terms in (2.22). □ 

The theorem may be applied to a variety of common discretizations. For example, 
consider an upwind discretization of the linear advection equation 

ut aux = 0, a > 0, 

on a grid with uniform spacing Aa;. In this example i = aui from which we observe 
that At2 = oo and Ati = ^Aa;. Assuming that the schemes forming the embedded 
pair have SSP coefficients C and C we find that the flux-based SPERK scheme is 
positivity-preserving provided 

At < -min(C,C)Ax. (2.24) 
a 

We can also apply the theorem to a discretization of the diffusion equation. Con- 
sider 

Ut - VUxx = 0, > 0, 

on a grid with uniform spacing Ax. Taking /i+i = ^(— Wi+i-t-Ui) gives Ati = At2 = 
^Ax^. Assuming that the schemes forming the embedded pair have SSP coefficients 
C and C, we find that the ffux-based SPERK scheme is positivity-preserving under 
the time step-size restriction 

At < ^min(C,C)Ax^ (2.25) 

The proof for equation-based partitioning follows in a similar manner. For com- 
pleteness we state the result below. 

Theorem 2.7. Suppose that the semi-discretization (2.1) of the one- dimensional 
conservation law (2.4) satisfies the positivity property 

W + AtG{W) > 
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for all W > and At < AtpE where the inequalities are taken component-wise. 
Further, suppose that an equation-based SPERK scheme is applied to (2.1), and that 
the schemes composing the embedded pair have SSP coefficients C and C. Then the 
full discretization is positivity-preserving for time steps satisfying 

0<At< mm{C,C)AtFE. 



3. Generalizations. In this section we provide some useful generalizations of 
our previous spatial partitioning methods. For notational simplicity we state these 
results for partitions that are invariant in time; however, it is also possible to vary 
the partitioning function x at each time level n. 

3.1. More than two methods. Both equation- and flux-based partitioning can 
be generalized to more than two methods. Accuracy and positivity results essentially 
identical to those in Sections 2.3 and 2.4 can be shown for both cases. 

3.1.1. Equation-based partitioning. Let {Ii,l2, ■ . .Tr) be a partitioning of 
the integers from 1 to iV (i.e., the partitioning satisfies Uj^.^j^Ifc = {1,2,..., N}, and 
Xj n Ifc = if 7^ fc). The generalization to r embedded methods defined by the 
coefficients A, . . . , b'--^^ is defined as 

s 

r s 

jjn+i ^ + AtJ2Yl diag(x^''^)G(r(^')), 

k=l = 1 

where 

(fc) _ f 1 if « e 

' 1 otherwise. 

3.1.2. Flux-based partitioning. Flux-based partitioning can also be written 
more generally in terms of fe^^^, . . . , 6^''). In this case we take 

/\t 

y(^) = [/"-_ ^ a,,i?<i>(r(^)), i<j<s, 
fc=i 

= [/" - i^diag(x('^))$(y«), 

to combine the r different embedded schemes. 

3.2. Blending and accuracy. Rather than shifting discontinuously between 
methods from one spatial point to the next, an appealing approach is to have a 
transition region in which a weighted average of the two methods is used, with the 
weight shifting from one method to the other over the transition region. The following 
proposition shows that this approach is still accurate in time, which is a necessary 
condition for the full SPERK discretization to be accurate. 

Proposition 3.1. Let {A, M^^), [A, fc^^)), ...,{A, 6'''') denote a set of embedded 
RK methods, and let pk denote the order of accuracy of method {A,b^^^). Then the 
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RK method 



r \ r 




has order of accuracy at least p = miiifc pk ■ 

Proof. The order conditions involve only expressions that are linear in the weights. 
Since all of the component method coefficients satisfy the conditions up to order 
p = minkPkj the averaged method does as well. □ 

Remark 3.1. The proposition applies only to embedded methods. If the compo- 
nent methods had different coefficient matrices A, their average would in general be 
only first-order accurate, even if they share common abscissas c. 

Theorems 2.1 and 2.4 for equation- and flux-based SPERK methods can be ex- 
tended to the blended case in this way. In practice, we observe the full accuracy when 
blending using either approach. 

4. Example: Spatially partitioned time-stepping for advection-diffusion. 

The semi-discretization of parabolic PDEs leads to systems of ODEs with eigenvalues 
near the negative real axis, whereas the semi-discretization of hyperbolic PDEs leads 
to systems with eigenvalues near the imaginary axis. Appropriate time integrators 
must include the corresponding portion of the complex plane in their absolute stabil- 
ity region. In problems with both hyperbolic and parabolic terms where the relative 
importance of these terms varies in space, e.g., in problems with strongly varying 
coefficients or mesh spacing, these criteria may be in conflict. 

Consider, for example, the linear, variable coefficient advection-diffusion equation 



with periodic boundary conditions on [0, 1] and initial conditions u{x, 0) = sin'^(27rx). 
Discretizing the spatial derivatives in (4.1) with three-point centered differences yields 
the flux-differencing method (2.5) with 



and a constant mesh spacing Ax. In our experiments, we take a{x) > and b{x) to 
be functions that are periodically defined by 

a{x) ^ + ^ arctan(1000(x - 0.1)) • arctan(1000(0.4 - x))) , -| < a; < |, 
b{x) = 1 + 38 (i + ^ arctan(1000(a; - 0.6)) • arctan(1000(0.9 - x))) , 3 < a; < |. 

As shown in Figure 4.1, there is a region centered at x ~ 0.25 where diffusion dom- 
inates and negative real axis inclusion is critical for the design of the time-stepping 
scheme. Similarly, in the region around x = 0.75, convection dominates and imag- 
inary axis inclusion is the most crucial design feature. The strong variation in the 
dominant term makes this a natural test problem for SPERK schemes. 

4.1. Second-order embedded pairs. We now design three- and four-stage 
embedded explicit Runge-Kutta pairs suitable for this problem. In each embedded 
pair, one method has a stability polynomial that maximizes (or nearly maximizes) the 
real axis interval of absolute stability, while the other method has a stability poly- 
nomial that maximizes (or nearly maximizes) the imaginary axis interval of absolute 



Ut -\- {b{x)u) 



X 



{a{x)u^) 



X 7 



(4.1) 
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Fig. 4.1. The variable coefficients a{x) (left) and b{x) (right). We observe a spatial variation 
in the relative importance of advection and diffusion. 



stability. Given a coefficient matrix the stability polynomial of a Runge-Kutta 
method (and hence its linear stability properties) can be chosen by solving a linear 
system of equations for the weights hj. We do not consider two-stage pairs since all 
second-order two-stage schemes have the same stability polynomial. 

The stability polynomial of a three-stage, second-order RK method has the form 
R{z) — l+z+z^ where as = Ac. For real-axis stability, we take the 3-stage 
Runge-Kutta-Chebyshev method of the family presented in [24], which we refer to 
as RKC(3,2), as the second scheme of our embedded pair. The stability polynomial 
of this method is a Bakker-Chebyshev polynomial and contains nearly the largest 
possible portion of the negative real axis over the class of three-stage second-order 
schemes. The first scheme is obtained by setting as = 1/4, which approximately 
maximizes imaginary axis inclusion, and selecting weights bj corresponding to the 
resulting polynomial. This gives the following embedded pair: 












3/8 


3/8 






3/8 


3/16 


3/16 




bT 


-1/3 


-20/9 


32/9 


bT 


-1/3 


4/9 


8/9 



The absolute stability regions are shown in Figure 4.2 (left). 

The stability polynomial of a four-stage, second-order RK method has the form 
R{z) = 1 -I- z -I- z^/2 -I- a^z^ + a^z^ where aa = Ac and — b^A^c. The optimal 
imaginary axis inclusion occurs for the choice as — 1/6, — 1/24, which gives the 
stability polynomial of the classical fourth-order Runge-Kutta method. Hence we 
take the classical method, which we shall denote by RK4, as the first scheme of our 
pair. Sufficient degrees of freedom remain to obtain a scheme with the same linear 
stability as the optimal four-stage, second-order scheme RKC(4,2). This gives the 
following embedded pair (with absolute stability regions as shown in Figure 4.2 right) 














1/2 


1/2 








1/2 





1/2 






1 








1 




bT 


1/6 


1/3 


1/3 


1/6 


bT 


2/125 


17/25 


36/125 


2/125 



(4.3) 



4.2. Numerical results. We now give the results of numerical experiments car- 
ried out using our three- and four- stage SPERK schemes. The partitioning parameter 
X is set equal to 1 wherever a{x) is more than 1/lOth of its maximum value; elsewhere 
it is zero. In all cases, we vary the time step-size and compute max-norm absolute 
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stability regions of 3-5taqe pair 




-6 -5 -4 -3 -2 -1 



Fig. 4.2. Absolute stability regions for the three-stage embedded pair (4.2) (left) and the four- 
stage embedded pair (4.3) (right). The shaded red region corresponds to the method optimized for real 
axis inclusion and the black line corresponds to the method optimized for imaginary axis inclusion. 




time step-size time step-size 



Fig. 4.3. Maximum norm absolute errors for SPERK schemes applied to the variable coeffi- 
cient advection- diffusion problem. The left plot gives results for three-stage schemes: (a) embedded 
scheme, (b) RKC(3,2), (c) equation-based partitioning, (d) flux-based partitioning. The right plot 
gives results for four-stage schemes: (a) RK4, (b) our variant of RKC(4,2) , (c) equation-based 
partitioning, (d) flux-based partitioning. 

errors by comparing to a highly accurate RK4 approximation of the spatially dis- 
cretized system at the final time t = 1. All computations use a constant mesh spacing 
of Ax = 1/200. 

The results for three-stage methods are given in the left-hand plot of Fig. 4.3. We 
find that RKC(3,2) is unstable over nearly the entire interval. The second method 
composing the embedded pair gives good results for At < 0.000125, but stability is 
lost for larger time step-sizes. The SPERK schemes derived from the combination are 
stable for values of Ai which are two or more times greater. Flux-based partitioning 
gives the smallest error of the three-stage schemes considered, and also allows for the 
largest stable time steps. 

The right-hand plot of Fig. 4.3 gives the results for four-stage methods. Similar to 
the three-stage case, our variant of RKC(4,2) is unstable for nearly the entire interval. 
The classical fourth-order Runge-Kutta method RK4 has the smallest error for At < 
0.00017 but becomes unstable for larger time step-sizes. Moving to SPERK schemes 
gives improved stability and second-order accurate results. Once again, we observe 
that flux-based partitioning has a smaller error and allows for larger stable time steps 
than equation-based partitioning. A review of the numerical results underlying the 
plot indicates that flux-based partitioning gives accurate results for At < 0.00065 and 
a regular, second-order error for At < 0.00031. 
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5. Example: Spatially partitioned time-stepping for WENO. Weighted 
essentially non-oscillatory (WENO) spatial discretizations involve an adaptive, data- 
dependent combination of several candidate stencils to compute fluxes [19]. The most 
commonly used scheme provides fifth-order accuracy in smooth regions of the solution 
and formally third-order spatial discretizations near shocks or other discontinuities. 
Assuming a positive flux function, the finite difference WENO scheme computes the 
positive flux differences in (2.5) with 

fj+h = '^o (i/j-2 - lfj-1 + ^fj) + wi (-g/j-^i + |/j + f/j+i) 

+w2(i/, + i/,+i-i/,+2), (5.1) 

where the weights ujq, wi, and uJ2 are chosen by computing data-dependent smoothness 
indicators [19] . We will use the WENO weights Ci;o,i,2 to select the mask x in a WENO- 
based SPERK scheme. 

5.1. An embedded pair for WENO. We start with the SSPRK(5,3) scheme 
[21]. This will be the lower-order scheme which is used in spatial regions where the 
solution is not smooth. If we embed this method in a larger 7-stage, 5th-order method, 
we have the following Butcher tableau 


















0.3773 


0.3773 












0.7545 


0.3773 


0.3773 










0.7290 


0.2430 


0.2430 


0.2430 








0.6992 


0.1536 


0.1536 


0.1536 


0.2385 






C6 


aei 


a62 


063 


064 


065 




C7 


an 


072 


073 


074 


075 


076 


fc^ 


0.2067 


0.2067 


0.1171 


0.1818 


0.2876 







bi 


b2 


63 


64 


65 


be h7 



(5.2) 



where we display only a few digits of the SSPRK(5,3) coefhcients. The unknown 
coefficients will determine an RK(7,5) scheme which will optimized for linear stability 
and use in spatial regions where the solution is smooth. 

A technique for determining these coefhcients by satisfying the order conditions 
is explored in [11] following a strategy designed in [23]. We attempt to maximize the 
linear stability properties of the RK(7,5) scheme; these are determined by the stability 
polynomial which in this case is parameterized by the coefficients of the and 
terms [6], polynomial expressions in the coefficients (5.2) which we denote by ag and 
ar respectively. We find that specifying the SSPRK(5,3) coefficients restricts the 
possible solutions to a straight line through the ag-a^ space. By examining the linear 
stability properties along this line (see Figure 5.1), we can choose one of the degrees 
of freedom (the "homogeneous polynomial" 1^5, see [11] for details) to maximize the 
linear stability properties of the resulting scheme. There are still six (nonlinear) order 
conditions to satisfy and six coefficients to be determined. Solving these remaining 
equations with a computer algebra system results in a single discrete solution and a 
one-parameter family of solutions. The one-parameter family was discarded because 
each member had either a negative node cg or unreasonably large coefficients. The lone 
solution did not have these deficiencies and we use it for our numerical experiments. 
The coefficients are shown to 15 digits in Table 5.1. 

5.2. Using the WENO weights to select a temporal scheme. We present 
some preliminary numerical results for the new RK(7,5) / SSPRK(5,3) embedded 
pair. The WENO weights wo,i,2 are used to determine which scheme is propagated 
for each grid point; if the weights are within some small threshold of their theoretical 
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-8.2 -0.15 -0.1-0.05 0.05 0.1 0.15 0.2 0.25 



Fig. 5.1. Choosing the linear stability properties for RK(7,5) scheme. Left: various mea- 
surements of the linear stability region versus a parameter /gs ■ Here p is the radius of the largest 
inscribed disc and p2 is the radius of the largest interval of the imaginary axis included. Right: 
the "WENO bean", the spectrum of the linearized WENO operator in smooth regions [8, 25, 13]: 
p3 measures the largest scaling of this bean that will fit in the linear stability region. We choose 
Iqs = 0.0037 and find pz 1.76, p ^ 1.58, and p2 ~ 1.2. 



021 


= 0.377268915331368, 










asi 


= 0.377268915331368, 


032 


= 0.377268915331368, 






04 1 


= 0.242995220537396, 


042 


= 0.242995220537396, 


043 


= 0.242995220537396, 


051 


= 0.153589067695126, 


052 


= 0.153589067695126, 


053 


= 0.153589067695126, 


054 


= 0.23845893284629, 


bi = 


= 0.206734020864804, 


62 = 


= 0.206734020864804, 


bs = 


= 0.117097251841844, 


64 = 


= 0.18180256012014, 


k = 


= 0.287632146308408, 


aei 


= 0.113015751552667, 


(162 


= 1.49947221487533, 


063 


= 0.134753400626063, 


064 


= -1.06421259296782, 


065 


= 0.205145170072233, 


071 


= -0.512110930783855, 


072 


= 3.91735780781337, 


073 


= -0.0470520461913835, 


074 


= -0.218621292015928, 


075 


= -1.64543995945252, 


0176 


= -0.494133579369683, 


bi = 


= 0.122097569374901, 


b2 = 


= 0.492898173466563, 


63 = 


= -0.232023614650883, 


bi . 


= -1.98394581022939, 


&5 = 


= 1.85394392181784, 


be = 


= 0.965538124667539, 


&7 = 


= -0.21850836444657. 



Table 5.1 

Coefficients of the embedded RK(7,5)/SSPRK(5,3) method. 



smooth-region values jq: and then we propagate the higher-order hnearly 
stable RK(7,5) solution. On the other hand, if the WENO weights do not indicate 
smoothness, then we propagate the lower-order SSPRK(5,3). The mask is widened 
slightly by four grid cells to increase the usage of the SSP scheme. We use a threshold 
of 0.06 and compute the mask based on [/" (the mask is fixed over each time-step). 

5.3. Numerical results. We perform convergence studies which demonstrate 
that the SPERK methods achieve the predicted orders of accuracy. Our tests use 
the inviscid Burgers' equation Ut -\- (u^)x = on the periodic domain [—1,1] with 
smooth initial conditions u{x, ^) ~ \ ^ \ cos (tt {x — ^^^^^^ ) ) ■ Table 5.2 shows that 
the methods achieve (at least) their expected orders of accuracy. Namely, five when 
the fifth-order scheme is used {x{x) = 1 for all x) and three when using the third- 
order scheme at any points. As noted in Section 3.2, x need not be a binary choice: 
we see that the method maintains third-order even when random values of x (in the 
range [0,1]) are selected at each point and at each time-step. Note also that the 
fiux-partitioned scheme exhibits an order of accuracy min(p,p) = 3, without the loss 
of one order as predicted by our Theorem 2.4; this is typical of what we observed in 
all our numerical tests. 

In addition to these finite difference convergence studies, we also tested the 
equation-based SPERK scheme with a finite-volume WENO code [12] and observed 
similar results. 
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error for different x(.^) functions 



partitioning 


Aa; 


1 









Heaviside(a;) 


random (0,1) 


equation- 


1/320 


1.29 X 10' 


-7 


5.59 X 10" 


-7 


1.29 X 10"^ 


3.33 X 10"^ 


based 


1/640 


4.09 X 10" 


-9 


7.00 X 10" 


-8 


4.44 X 10"^ 


3.76 X 10~* 




1/1280 


1.29 X 10" 


'10 


8.64 X 10" 


-9 


5.48 X W-^° 


4.39 X lO"'' 


est. order 




4.99 




3.02 




3.45 


3.08 


flux-based 


1/320 


2.01 X 10" 


-8 


1.08 X 10" 


-7 


2.44 X 10"* 


6.63 X 10"** 




1/640 


6.38 X 10" 


-10 


1.34 X 10" 


-8 


1.87 X 10"^ 


8.06 X lO"'* 




1/1280 


2.00 X 10" 


-11 


1.65 X 10" 


-9 


2.18 X 10"^" 


1.05 X 10"^ 


est. order 




4.99 




3.02 




3.44 


2.99 



Table 5.2 

Convergence studies for various choices of the mask function x{^)t performed on Burgers' 
equation with smooth initial conditions discretized using WENO and the SPERK 5th- order /3rd- 
order pair. We use a CFL number of 1.2 and the error is measured in an approximate L2 norm. 



method 




time 


error 


error {Ax/2) 


est. order 


TV increase 


RK(7,5) 


t 


= 0.25 


2.01e-8 


6.38e-10 


4.98 







t 


= 1.25 


8.72e-3 


8.79e-3 




0.264 


SSPRK(5,3) 


t 


= 0.25 


1.08e-7 


1.34e-8 


3.01 







t 


= 1.25 


1.35e-3 


9.53e-4 


0.50 





SPERK 


t 


= 0.25 


2.01e-8 


6.38e-10 


4.98 







t 


= 1.25 


1.35e-3 


9.53e-4 


0.50 






Table 5.3 



The flux-partitioned SPERK scheme combined viith WENO has the best features of each method. 
Here the initial smooth curve sharpens into a shock. At small times t = 0.25 before shock formation, 
both of the underlying schemes exhibit their design orders of 5 and 3 respectively. At a later time 
t = 1.25, a shock has formed; all schemes exhibit larger errors as the solution is no longer smooth. 
However, RK(7,5) now exhibits spurious oscillations, indicated by an increase in the total variation 
seminorm. We use a CFL number of 1.2, Ax' = gig, and error measured in an approx. L2 norm. 



5.3.1. Total variation tests. We consider the discrete total variation seminorm 
of the solution on a test problem consisting of the inviscid Burgers' equation with a 
square wave initial condition in [0,1]. We determined the largest CFL a, At = 
aAx such that the resulting solution experiences no increase in the TV seminorm. 
For the flux-partitioned scheme we found that the SSPRK(5,3) scheme (i.e., x = 
everywhere) is no longer TVD for a > 1.3. The RK(7,5) solution is no longer TVD for 
cr > 0.5. With the SPERK embedded pair, the loss of the TVD property occurs for 
a > 1.3 (that is, the same as the SSPRK(5,3) scheme). The results using equation- 
based partitioning are about the same. 

The example in Table 5.3 demonstrates that the RK(7,5)/SSPRK(5,3) SPERK 
scheme offers high-order accuracy when the solution is smooth and remains TVD 
when the solution is non-smooth; the best of both worlds! 

6. Discussion and Future Work. This paper considers spatially partitioned 
embedded Runge-Kutta (SPERK) schemes. Such methods give an efficient means 
to apply two or more Runge-Kutta methods to a spatially discretized PDE. In this 
paper, SPERK schemes are applied in two ways. The first of these partitions the 
ODEs, and is therefore referred to as equation-based partitioning. In these meth- 
ods the order of accuracy is shown to equal the minimum of the order of accuracy 
of the schemes composing the embedded pair. Equation-based partitioning is not 
conservative when applied to a conservation law, however. The second partitioning 
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method partitions by fluxes. This flux-based partitioning has the advantage of being 
conservative when applied to a conservation law. Theoretically, it may lead to a loss 
of one order of accuracy when compared to the underlying schemes; however, this 
loss of accuracy is not observed in practice. We show that both equation- and flux- 
based partitioning are positivity-preserving under a suitable time step restriction when 
the underlying schemes are strong-stability-preserving (SSP). Numerical experiments 
on a spatially discretized variable coefficient advection-diffusion equation show that 
SPERK schemes can have superior linear stability characteristics than either scheme 
individually. SPERK schemes may be applied to weighted non-oscillatory (WENO) 
spatial discretizations of conservation laws. Our partitioning is carried out using the 
weights introduced in the WENO spatial discretization. Specifically, a fifth-order 
Runge-Kutta method is used in smooth regions where WENO chooses a fifth-order 
stencil, and a third-order SSP Runge-Kutta method is used in non-smooth regions 
where WENO chooses a formally third-order stencil. We find that this combination 
avoids oscillations in problems with shocks, and gives fifth-order accuracy in smooth 
flows. 

As part of this paper, several explicit SPERK schemes are designed and numeri- 
cally validated. Future work will carry out a more systematic development of SPERK 
schemes. Of particular interest for us are methods that have good monotonicity where 
the solution is nonsmooth, and good accuracy elsewhere. Such methods are partic- 
ularly attractive for approximating conservation laws with isolated shocks and other 
nonsmooth features. 

Schemes with improved local absolute stability properties (like those designed in 
Section 4) seem promising for problems with strong spatial variation in the coefficients 
or grid size. The stability polynomial optimization algorithm developed in [10] could 
be used more generally to design appropriate methods for such problems. 

Implicit SPERK schemes are also of interest. For example, in regions with nons- 
mooth features, one may wish to use backward Euler since it is unconditionally mono- 
tone, and therefore avoids spurious oscillations. However, in regions characterized by 
smooth solutions, a higher-order method may be preferred since backward Euler is 
only first-order accurate and is highly dissipative. The corresponding SPERK com- 
bination may give improved accuracy and less smearing than backward Euler, while 
simultaneously offering improved monotonicity over second- or higher- order schemes. 
A detailed examination of implicit, equation-based partitioning based on combina- 
tions of the ^-scheme appears in the report of Van Slingerland, Borsboom and Vuik 
[15]. They find optimal local 0-values to minimize the numerical diffusion, while 
ensuring that the scheme is stable, positivity-preserving, and non-oscillatory. It is 
natural to consider extensions of their work and our own. For example, the develop- 
ment and analysis of higher-order implicit SPERK schemes as well as the study of 
implicit methods based on flux partitioning would be particularly interesting. 

Theoretically, the accuracy of flux-based partitioning is expected to be less than 
the accuracy of the schemes composing the embedded pair but this loss of accuracy 
is not seen in our tests. This effect may be similar to what is observed in multirate 
schemes. For example, the Osher-Sanders scheme [14] is locally inconsistent, but has 
been found to give flrst-order convergence in practice. An explanation for this is given 
in [7] where it is shown that local inconsistencies may not show up in the global errors 
due to cancellation and damping effects. We hope to transfer this analysis to the 
case of SPERK schemes to obtain a better understanding of the unexpectedly good 
performance of flux-based partitioning. 
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